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Abstract. Data assimilation procedures have been developed for thermospheric mod¬ 
els using satellite density measurements as part of the EU Framework Package 7 AT- 
MOP Project. Two models were studied; one a general circulation model, TIEGCM, and 
the other a semi-empirical drag temperature model, DTM. Results of runs using data 
assimilation with these models were compared with independent density observations from 
CHAMP and GRACE satellites throughout solar cycles 23 and 24. Time periods of 60 
days were examined at solar minimum and maximum, including the 2003 Hallowe’en storms. 
The differences between the physical and the semi-empirical models have been charac¬ 
terised. Results indicate that both models tend to show similar behaviour; underestimat¬ 
ing densities at solar maximum, and overestimating them at solar minimum. DTM per¬ 
formed better at solar minimum, with both models less accurate at solar maximum. A 
mean improvement of ~ 4% was found using data assimilation with TIEGCM. With fur¬ 
ther improvements, the use of general circulation models in operational space weather 
forecasting (in addition to empirical methods currently used) is plausible. Future work 
will allow near-real-time assimilation of thermospheric data for improved forecasting. 


1. Introduction 

The Advanced Thermosphere Modelling of Orbital Pre¬ 
diction (ATMOP) project was designed to provide a Euro¬ 
pean capability for nowcasting and forecasting of the ther¬ 
mosphere. Changes in thermospheric density affect the drag 
experienced by low-orbiting satellites and hence their orbits. 
The resulting loss in orbit predictability is problematic for 
space agencies and satellite operators, with thousands of 
objects orbiting the Earth and the majority of them in low- 
earth orbit [see review by Vallado and Finkleman, 2014]. 
Space weather is one cause for these density changes, as 
the associated changes in solar or geomagnetic activity af¬ 
fects the coupled thermosphere-ionosphere system [Rees and 
Fuller-Rowell, 1989; Ftuba et ai, 2014]. 

A major focus of ATMOP was to investigate the ben¬ 
efit of data assimilation to predictive models. While this 
is a standard method in meteorological forecasting, it is a 
relatively new approach within the space weather commu¬ 
nity. This paper compares both semi-empirical and physical 
models using this new approach. Semi-empirical models are 
more rapid to run, and typically outperform physical mod¬ 
els, hence are currently used for operational orbit computa¬ 
tions. However, the physical models may perform better 
in unusual conditions, such as those encountered in geo¬ 
magnetic storms. Two models are the focus of this work 
- the first a version of the semi-empirical Drag Tempera¬ 
ture Model ]DTM-2012; Bruinsma et ai, 2003], and the sec¬ 
ond the physical Thermosphere Ionosphere Electrodynam¬ 
ics General Circulation Model [TIEGCM; Richmond et ai, 
1992]. 

A large amount of previous work has examined the rel¬ 
ative accuracy of thermospheric general circulation models 
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compared to empirical methods. For example, Anderson 
et al. [1998] intercompared five models including TIEGCM, 
in order to determine why several physical models consis¬ 
tently underestimate the ionospheric F region peak electron 
density. Qian and Solomon [2012] examined temporal and 
spatial variations of thermospheric density via a compari¬ 
son of TIEGCM with an empirical model, NRLMSISE-00 
[Picone et al, 2002]. In-situ observations from satellites are 
used in these comparisons to determine whether the mod¬ 
els accurately represent thermospheric conditions [see also, 
Buonsanto et al., 1997; Pavlov et al, 2000; Forbes et al., 
2005; Sutton et al., 2005; Tsagouri et al., 2013]. Readily 
available satellite data are vital for space weather opera¬ 
tions, being used in data assimilation techniques to improve 
the thermospheric models that are used to predict the drag 
experienced by low-orbiting spacecraft [Fuller-Rowell et al., 
2006]. 

A typical atmospheric model set up with boundary con¬ 
ditions closely matching reality will nevertheless produce 
predictions which gradually diverge from reality due to in¬ 
complete characterisation of the physics and chaotic effects. 
Thermospheric models are also strongly influenced by exter¬ 
nal drivers; geomagnetic and solar drivers (which are gen¬ 
erally represented by way of input indices to the models), 
and the lower and upper boundaries implemented to repre¬ 
sent the interaction with the lower atmosphere below and 
plasmasphere above. By assimilating observations into a 
physical model at a given time, model fields at that time 
can be brought closer to the real data values. Bringing the 
model state closer to reality reduces the divergence of the 
ensuing prediction from reality, thereby increasing the accu¬ 
racy of the model forecast. Data assimilation repeats this 
process regularly, to provide a running forecast with better 
overall accuracy. 

Much previous thermospheric data assimilation work has 
focused on using Kalman filters [Fuller-Rowell et al., 2004], 
particularly the development of ensemble Kalman filtering 
methods ]EnKF; Codrescu et al., 2004; Matsuo et al., 2012] 
with recent advances in computer technology. Matsuo et al. 
]2012] found the root-mean-square (RMS) error of model- 
predicted neutral mass density could be reduced by up to 


1 



X - 2 


MURRAY ET AL.: THERMOSPHERIC MODELLING DURING SOLAR CYCLES 23 AND 24. 


50% in the case that as many as nine ensemble members 
are used. Codrescu et al. [2004] found the error of their 
model initial state reduced from ~25% to ^10% using an 
EnKF technique with ten members. Matsuo et al. [2013] 
used EnKE with inferred neutral density satellite observa¬ 
tions to improve a model neutral density specification in the 
vicinity of the satellite. Matsuo et al. found that a global 
impact was achieved if accompanied by the estimation of the 
primary driver of the density variable (such as solar EUV 
flux), highlighting the importance of external forcing to the 
system. 

Two data assimilation methods were developed in AT- 
MOP; one using DTM-2012, and the other TIEGCM. A ver¬ 
sion of DTM was developed that can assimilate total density 
data in near-real time, which will make orbit predictions sig¬ 
nificantly more accurate. Data assimilation techniques were 
also developed for use with TIEGCM with the aim to create 
a more physically accurate global analysis and forecast sys¬ 
tem for the thermosphere. While the assimilation methods 
developed are described in detail elsewhere [see Bruinsma 
et al, 2012, and Henley et al, manuscript in preparation], 
this paper provides first results of an intercomparison of the 
newly developed models. Both TIEGCM and DTM-2012 
are compared with observations taken at various periods 
throughout solar cycles 23 and 24 in order to determine if a 
more complex physical model with data assimilation could 
be more accurate for forecasting efforts than the currently- 
used empirical methods. This validation effort is crucial 
to determine if these new techniques are beneficial for cur¬ 
rent space weather operations. The observations and models 
used for the comparison will be outlined in Section 2, and 
methods used for analysis described in Section 3. The re¬ 
sults of the study will be presented in Section 4, while some 
conclusions of the findings and possible future work will be 
discussed in Section 5. 

2. Observations and Models 

Atmospheric drag on satellites varies strongly as a func¬ 
tion of thermospheric mass density [for further discussion 
see Qtan and Solomon, 2012; Vallado and Finkleman, 2014]. 
This work uses in-situ observations of atmospheric densities 
in the upper thermosphere for both the data assimilation 
techniques and the validation and intercomparison of the 
model results. The thermospheric densities were inferred 
from accelerometers that fly on low-earth-orbiting satellites, 
namely the CHAllenging Minisatellite Payload [CHAMP; 
Wickert et al., 2001] and the Gravity Recovery And Climate 
Experiment [GRACE; Tapley et al., 2004]. See Bruinsma 
et al. [2004] for specific details on how the densities are cal¬ 
culated. The height at which the densities are available 
depends on the satellite orbital height; this is ^450-300 km 
for CHAMP and ~490-410 km for GRACE (as of 2014 Oc¬ 
tober). The height value is dependant on the time period of 
observations; CHAMP data are available from 2000 to 2010, 
and GRACE data from 2002 to present. The satellites ob¬ 
serve at ^15 orbits per day at all 360° of longitude, with a 
spacing of 25° at the equator per day, and a latitude reso¬ 
lution of 0.5-1°. The horizontal resolution is 80 and 40 km 
along the orbit for CHAMP and GRACE, respectively. Ob¬ 
servations are made every 10 s with CHAMP and 5 s with 
GRACE. The satellite data have a mean and RMS observa¬ 
tion error both ranging between 1-10%, which is a function 
of latitude and geomagnetic activity. 

The drag temperature model used for comparison in this 
paper, DTM-2012 (hereafter DTM), is an improved ver¬ 
sion of DTM-2009 [Bruinsma et al., 2012] that is trained 
on a larger data set. It has been constructed using the 
full CHAMP and GRACE high-resolution accelerometer in¬ 
ferred density data sets, besides new mean total density data 
as well as historical mass spectrometer data. The 81-day 


mean and 1-day delayed solar radio flux at 10.7cm (F10.7) 
are used as solar inputs, and the Am geomagnetic index is 
also used. DTM is constructed by fitting to the underlying 
density database, as good as possible in the least-squares 
sense, to reproduce the mean climatology of the thermo¬ 
sphere [see Bruinsma et al, 2004]. The result includes a 
density value at a particular latitude and longitude between 
200-1000 km, with up to 1° resolution. Errors on short time 
scales of a few days or less can be of the order of tens of per¬ 
cent. A more detailed description of the updated DTM with 
data assimilation can be found in Bruinsma et al. ]2012]. 

The physical model used in this work, TIEGCM, is a first- 
principles, three-dimensional, non-linear representation of 
the coupled thermosphere and ionosphere system. TIEGCM 
uses spherical geographic coordinates, with latitude ranging 
from —87.5° to 87.5° in 5° increments, longitude ranging 
from —180° to 180° in 5° increments, and a 60-second time 
step. The lower boundary at -^97 km extends up to -^500- 
700 km depending on solar activity. The migrating and 
semi-diurnal tides are specified at the lower boundary using 
the Global Scale Wave Model [Hagan et al., 1995], which 
does not consider the effects of planetary waves and non¬ 
migrating tides. The vertical coordinate is a log-pressure 
scale, ln(po/p), where p is pressure and po is a reference 
pressure, and pressure levels range from —7 to 7 increasing 
in half-scale-height increments. 

Neutral density observations were assimilated into the 
model using an ensemble optimal interpolation method 
]EnOI; Oke et al, 2002; Evensen, 2003]. EnOI uses an 
ensemble approach to help determine how much trust to 
put in the model prediction at any given time (as opposed 
to how much trust to put in the observations). An ensem¬ 
ble of nine models are run offline, with each model regularly 
perturbed with smoothed random temperature fluctuations. 
A background density error covariance matrix is then cal¬ 
culated from these non-assimilative, independent ensemble 
members. This is combined with the background and obser¬ 
vations to produce a density analysis. The analysis is then 
converted to a temperature analysis (as density is a derived 
field in TIEGGM), which is used to initialise the model run. 
This is the first time an EnOI method has been used for 
thermospheric data assimilation, with much previous work 
focusing on EnKF (see Section 1). For further details on 
the assimilation scheme, see Henley et al. (manuscript in 
preparation, 2015). 

For consistency with DTM input parameters, TIEGCM 
was run using the Heelis high-latitude ion convection model 
[Heelis et al., 1982] as a magnetospheric input, which uses 
F10.7 solar flux and Kp planetary geomagnetic index as in¬ 
puts. The model provides various advantages for develop¬ 
ment, speed, robustness, and scientific return compared to 
empirical modelling. Results from TIEGCM have been pre¬ 
viously compared favourably with CHAMP data on a large 
scale [Sutton, 2008; Qian and Solomon, 2012], although it 
has also been noted that TIEGCM RMS errors typically 
gradually increase with a decline in solar activity [Kim, 
2011 ]. 

3. Comparison 

TIEGCM and DTM results were compared with CHAMP 
and GRACE observations at three periods throughout so¬ 
lar cycles 23 and 24. Solar equinox periods were selected 
for ease of comparison of model-run results. For solar min¬ 
imum, a 60-day period was chosen beginning 2009 March 
15, and for solar maximum another 60-day period begin¬ 
ning 2003 March 15 was selected. To examine the impact of 
more severe stormy conditions, a final 60-day period was ex¬ 
amined beginning 2003 October 15, which encompasses the 
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Figure 1. DTM (upper row) and TIEGCM (middle row) 
density values on 2009 March 01 00:00UT, and CHAMP 
(lower row) density values between 2009 March 01 00:00- 
01:30UT. The DTM and CHAMP values have been in¬ 
terpolated to TIEGCM pressure level 21 at ~ 10“® Pa 
(~335km). 


Hallowe’en storms. These storms are defined as a sequence 
of events between 2003 October 28 to November 4, during 
which massive solar flares and geomagnetic storms occurred. 
The impact of this event on the near-Earth environment has 
been heavily studied, including analysis of observations [e.g., 
Thomson et al, 2004; Mannucci et al, 2005; Bruinsma et al., 
2006] as well as the use of simulations [e.g., Toth et al., 2007; 
Manchester et al., 2008]. Gopalswamy et al. [2005] gives a 
detailed review of the events, while Oler [2004] outlines the 
prediction performance of space weather forecasting centres 
worldwide following these events. 

A large database of observations were assimilated into 
DTM (see Section 2) before running the model on the time 
periods selected. TIEGCM was run with 1-hourly data as¬ 
similation using the EnOI system (see Section 2) along with 
CHAMP observations for these periods. It is worth not¬ 
ing the impact of assimilation is at a local level since the 
satellite data give density values at a particular latitude, 
longitude, and height, and a localisation scheme is also ap¬ 
plied. The CHAMP data are decimated from their original 
10 s cadence to 120 s cadences in order to help stability, 
and this also provides semi-independent CHAMP data sets 
to use for validation. More formal validation is provided by 
independant GRACE observations. 

The model results were interpolated to the altitude, lati¬ 
tude and longitude of the particular spacecraft at that time 
for comparison purposes. The interpolation procedure used 
was the same as that used within the TIEGCM data as¬ 
similation code for consistency (Henley et ah, manuscript 


in preparation, 2015). In order to accurately compare the 
output of the two models, DTM was run for the particular 
altitude, latitude, and longitude of the spacecraft at time of 
observation. DTM densities were calculated every 5 or 10 
seconds (depending whether CHAMP or GRACE was being 
compared), while TIEGCM resulting densities were saved 
every 15 minutes. It is worth noting that both CHAMP 
and GRACE spacecraft have different local times for the pe¬ 
riods studied here, the difference ranging between ~l-2 LT. 
It useful to compare the model results to these observations 
for validation purposes since they have an altitude differ¬ 
ence of ~100-150km for the time periods studied. However, 
a more in-depth comparison of latitudinal or longitudinal 
variation is not undertaken here. 

Typical model outputs are presented in Figure 1. A 2D 
map of DTM density at 2009 March 01 OOiOOUT is shown in 
the upper row, with an equivalent TIEGCM map in the mid¬ 
dle row. CHAMP densities are also plotted in the lower row 
for an entire 90-minute orbit. Note the CHAMP spacecraft 
was at an altitude of ~328km during this orbit, however 
the DTM and CHAMP densities have been interpolated to 
TIEGCM pressure level 21 (1.18 x 10“® Pa) for ease of com¬ 
parison. TIEGCM tends to have a narrower range of density 
values compared to DTM. It is also worth noting that both 
models overestimate the density in the points corresponding 
to the CHAMP measurements. This is a result that is emu¬ 
lated throughout solar minimum, as will be discussed in the 
following section. 

4. Results 

The results of the thermospheric model comparison for 
solar minimum are presented in Figure 2. The lower mid¬ 
dle row shows the F10.7 solar flux (magenta) and Kp indices 
(black) that were used as inputs to both models. With F10.7 
virtually constant at ~70, and Kp generally below ~3, this 
can be considered a quiet period in terms of storms, and 
thus a good indicator of results at typical solar minimum 
conditions. In soft X-ray range, solar flares are classified as 
A-, B-, C-, M- or X- class according to peak flux measured 
near Earth by the Geostationary Operational Environmen¬ 
tal Satellite (GOES) over 1-8 A in Watts m Solar ac¬ 
tivity was extremely low throughout this period, with no 
solar flares greater than B-class, as indicated by the GOES 
flux plot in the lower row of Figure 2. Two B-class flares 
occurred on 2009 March 26 (see lower row of Figure 2), with 
no major active regions on the solar disk at the time. A 
larger number of A- and B- class flares occurred at the end 
of the period, again with no major active regions on the 
solar disk. 

The upper row of Figure 2 shows the orbit-averaged den¬ 
sity values of GHAMP (black), TIEGGM (blue), and DTM 
(green). The upper middle row shows the density difference 
between GHAMP and TIEGCM (blue), and CHAMP and 
DTM (green). The results indicate that DTM performs bet¬ 
ter at solar minimum, with its resulting densities closest to 
the actual observations from CHAMP. However, TIEGCM 
also performs very well, with densities only slightly higher 
than DTM overall. The density values are relatively sta¬ 
ble throughout the period. Both models mirror any gradual 
changes in density, with DTM mirroring rapid changes par¬ 
ticularly well during solar minimum conditions. 

The upper row of Figure 5 shows the same results as Fig¬ 
ure 2, however for a shorter time period of 12 hours on 2009 
March 15. This zoom-in reveals greater detail in the evolu¬ 
tion of the densities, showing typical motions of the CHAMP 
spacecraft as it changes in altitude, latitude, and longitude. 
Both models mirror the CHAMP trend well, with DTM per¬ 
forming better than TIEGCM as the spacecraft crosses the 
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Figure 2. Upper row: Orbit-averaged CHAMP (black), DTM (green), and TIEGCM (blue) densities 
for 2009 March 15 to 2009 May 13. Upper middle: Density difference between CHAMP observations and 
the two models. Lower middle row: F10.7 (magenta) and Kp (black) values used as inputs to the models. 
Lower row: GOES peak flux (red), where the dashed horizontal lines intersecting the right y-axis indicate 
flare class. The altitude of the CHAMP spacecraft during this period varied between ~350-316 km. 


night-side. It is clear from both figures that the models gen¬ 
erally overestimate the density during solar minimum, with 
resulting model densities consistently higher than observa¬ 
tions throughout the period of analysis. 

The same analysis was performed on CHAMP and both 
models at solar maximum, the results of which are shown in 
Eigure 3. It is clear from the lower middle row of the figure 
that there were more active conditions during this period, 
with F10.7 solar flux varying between ~100-160, and Kp 
index reaching ~7. Although no significantly large storms 
occurred during this time, the values are sufficiently high to 
contrast with the solar minimum results. Solar activity was 
much higher, with consistent flaring throughout the period, 
including multiple large flares (see lower row of Figure 3). 
For example, a number of M- and X- class flares occurred be¬ 
tween 2003 March 17-19 from a very complex active region 
just west of solar disk centre. This region was classified as 
according to the modified Mount Wilson sunspot classi¬ 
fication scheme [Bray and Loughhead, 1964]. The impact of 
these eruptions is not noticeable from the F10.7 and Kp val¬ 
ues at this time, with no significant Earth-directed coronal 
mass ejections associated with the flares. 

As indicated by the orbit-averaged densities and differ¬ 
ence shown in the upper and upper middle rows of Eigure 3, 
TIEGCM generally outperforms DTM at solar maximum. 
DTM does mirror increases in the CHAMP densities better 


than the physical model at the sharp increase seen on 2003 
March 20 in CHAMP density measurements (likely related 
to the increased solar activity at this time). DTM also shows 
a sharp increase on April 30 that mirrors a smaller increase 
in CHAMP. TIEGCM densities do not change significantly 
on March 20, but do slightly increase on April 30. How¬ 
ever, DTM also shows sharp increases/decreases in values 
during periods of increased Kp index that are not reflected 
in CHAMP or TIEGCM densities (e.g., April 16 and May 
9). This suggests DTM is more sensitive to sharp changes 
in input parameter values than TIEGCM is. The upper row 
of Figure 5 shows the solar maximum results over a shorter 
12-hour period on 2003 March 15. This confirms that both 
models generally underestimate thermospheric density at so¬ 
lar maximum (in contrast to solar minimum results), and 
TIEGCM generally performs better during this period. 

The results of analysis during the 2003 October and 
November storms are shown in Figure 4. The signifi¬ 
cantly stormy conditions that occur towards the end of 
2003 October are evident from GOES soft X-ray flux 
in the lower row, with multiple X-class flares, including 
an X17.2 flare on October 28. This resulted in a G5 
level (http://www.swpc.noaa.gov/NOAAscales/) geomag¬ 
netic storm at Earth that lasted two days. The source 
regions of these large eruptions are discussed in detail by 
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Figure 3. Upper row: Orbit-averaged CHAMP (black), DTM (green), and TIEGCM (blue) densities 
for 2003 March 15 to 2003 May 13. Upper middle: Density difference between CHAMP observations and 
the two models. Lower middle row: F10.7 (magenta) and Kp (black) values used as inputs to the models. 
Lower row: GOES peak flux (red), where the dashed horizontal lines intersecting the right y-axis indicate 
flare class. The altitude of the CHAMP spacecraft during this period varied between ^433-399 km. 


Zhang et al. [2003]. It is worth noting that the largest flare 
ever measured occurred on November 4, believed to be as 
large as X45 [Thomson et al, 2004). However, this event 
was not Earth-oriented and only resulted in high-latitude 
auroras. The sequence of events are reflected in the large 
Kp and F10.7 values shown in the lower middle row of Fig¬ 
ure 4. A second large peak in Kp occurs on November 20, 
however only M-class flares occurred during this period and 
El0.7 did not increase as much as previously. 

The density values and differences shown in the upper and 
upper middle rows of Figure 4 reflect the particularly active 
conditions, with two large peaks in density at the end of Oc¬ 
tober and November. Both models perform well in mirroring 
the observed increased densities, however they both under¬ 
estimate the values again, particularly at the large peaks. 
TIEGCM more accurately follows the changes in CHAMP 
values, with DTM showing a more sporadic variance. To¬ 
wards the end of the time period DTM values sharply rise, 
which is not reflected in the observed values. However, this 
rise does occur soon after an increase in Kp. The zoomed-in 
plots in the lower row of Figure 5 confirm TIEGCM be¬ 
ing the more accurate model in this case, since it generally 
follows changes in CHAMP values better than DTM. Com¬ 
pared to results from 2003 March, it seems that TIEGCM 
reacts best to sharp rises in Kp values greater than ~ 8, and 
does not react as significantly to changes in more moder¬ 


ate values. Carter et al. [2014] found the most important 
source of variability in TIEGCM originates from the Kp in¬ 
dex, so it is unsurprising that it the model results reflect 
sharp changes in this parameter. 

GRACE observations were also available for two time pe¬ 
riods at solar minimum (2009 March) and maximum (2003 
October). The same analysis as above was performed us¬ 
ing GRACE data for these periods, as well as using the 
data for further analysis of the previous results. The re¬ 
sulting densities from running TIEGCM with assimilation 
of CHAMP data were interpolated to the altitude, latitude, 
and longitude of GRACE at the time of the observation in 
order to gain a better insight into the overall impact of run¬ 
ning TIEGCM with assimilation using such a small amount 
of data. This ensures no bias with interpolating CHAMP- 
assimilated TIEGCM results to CHAMP observations (and 
then comparing with same). 

Results from the solar minimum period beginning 2009 
March 15 (upper and upper middle rows) and the stormy 
period beginning 2003 October 15 (lower middle and lower 
rows) are shown in Figure 6, including DTM runs for com¬ 
parison. DTM was run for the particular time, altitude, 
latitude, and longitude of GRACE as previously. Com¬ 
paring the results for the solar minimum period in Fig¬ 
ure 2, TIEGCM performs similarly within the same er¬ 
ror range, however DTM more significantly outperforms 
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Figure 4. Upper row: Orbit-averaged CHAMP (black), DTM (green), and TIEGCM (blue) densities 
for 2003 Ocotber 15 to 2003 December 13. Upper middle: Density difference between CHAMP observa¬ 
tions and the two models. Lower middle row: F10.7 (magenta) and Kp (black) values used as inputs to 
the models. Lower row: GOES peak flux (red), where the dashed horizontal lines intersecting the right 
y-axis indicate flare class. The altitude of the CHAMP spacecraft during this period varied between 
~422-386 km. 


TIEGCM compared to previously. A discrepancy occurs on 
2009 May 3 when both models show a slight increase com¬ 
pared to a drop in observed density values, however this is 
likely due to an issue with the GRACE measurements for 
that particular date rather than the model results. It is also 
interesting to note that the models underestimate the ob¬ 
served GRACE values during this period, while they overes¬ 
timated CHAMP values previously. Comparing the results 
from the Hallowe’en storm period in Figure 6 to those shown 
in Figure 4, this time TIEGCM outperforms DTM, as it did 
when compared to CHAMP observations. TIEGCM density 
values match GRACE observations closer than DTM, which 
has a large number of outliers during the stormy peaks at 
the end of 2003 October and 2003 November. 

All results from the model comparisons are summarised 
in Table 1, showing monthly mean percentage difference be¬ 
tween measured satellite densities and the resulting densi¬ 
ties, from various model runs. As found above, the table 
confirms that DTM performs best compared to CHAMP 
at solar minimum, followed closely by TIEGCM, with both 
models slightly overestimating the CHAMP density val¬ 
ues. DTM densities are much closer to GRACE observa¬ 
tions compared to TIEGCM during this period, with both 
models underestimating values at the higher satellite alti- 


Table 1. Mean percentage difference between measured 
CHAMP (upper table) and GRACE (lower table) densities 
and resulting densities from various model runs. Values were 
calculated for the time periods beginning in March 2003, Oc¬ 
tober 2003, and March 2009. 


Validating 

Data 

Model Run 

Mean % Difference 

March ‘03 

Oct. ‘03 

March ‘09 

CHAMP 

TIEGCM 

30.09 

9.76 

-26.86 


TIEGCM (DA) 

25.50 

5.58 

-23.35 


DTM 

51.99 

19.27 

-21.92 

GRACE 

TIEGCM 

- 

10.66 

46.66 


TIEGCM (DA) 

- 

6.60 

43.58 


DTM 

- 

43.33 

5.54 


tude. There is a clear difference during solar maximum, 
with TIEGCM being more accurate at these times, and 
both models underestimating the CHAMP density associ¬ 
ated with higher activity. The sharp changes in density dur¬ 
ing the Hallowe’en storms result in values being smaller dur¬ 
ing this period, however TIEGCM again outperforms DTM 
here, most significantly when compared with the GRACE 
observations (as is clear from Figure 4, DTM overestimates 
the density values at the peaks). 
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Figure 5. Left column: CHAMP, DTM, and TIEGCM densities over 12 hours. Right column: Density 
difference between CHAMP observations and the two models. The upper to lower rows show results 
for 12 hours on 2009 March 15, 2003 March 15, and 2003 October 15, respectively. The altitude of the 
CHAMP spacecraft during these periods varied between ~346-317 km, ~432-401 km, and ~422-391 km. 


The table also notes the differences between running 
TIEGCM with and without data assimilation, as it is worth 
determining how much of an improvement the assimilation 
system made to the physical model. Similar comparisons for 
DTM are presented in Bruinsma et al. [2012]. The results 
indicate that using the EnOI data assimilation technique has 
a small positive impact on the TIEGCM results, improving 
performance by ~4% overall when a 1-hour assimilation cy¬ 
cle is used. The data assimilation technique improves the 
results greater during periods of solar maximum than mini¬ 
mum, the largest improvement found during the Hallowe’en 
storms. 

5. Discussion and Conclusions 

This paper has characterised the relative merit of two new 
modelling approaches with data assimilation capabilities de¬ 
veloped during the ATMOP project; the physical model, 
TIEGCM, and the semi-empirical model, DTM. This has 
been done for 60-day periods for equinox at solar minimum 
(from 2009 March), solar maximum (from 2003 March), and 
the Hallowe’en 2003 storms (from 2003 October). Model 
results were validated against satellite data to investigate 
whether the improved physical model can outperform the 
semi-empirical model that is currently used operationally. 


The solar minimum equinox results show that the mod¬ 
els represent the thermospheric density well, with the DTM 
model being closest to CHAMP and GRACE observations, 
followed closely by TIEGCM. Interestingly, both models 
typically overestimate the CHAMP density, but underes¬ 
timate the GRACE density. It is worth noting that satel¬ 
lite drag data studied by Solomon et al. [2011] indicates 
that the thermosphere was lower in density, and therefore 
cooler, during the protracted solar minimum period of 2007- 
2009 than at any other time in the past 47 years. Thus 
the overestimated values compared to CHAMP here may be 
representative of this particular ‘low’ rather than represen¬ 
tative of general solar minimum conditions - other periods 
would need to be examined for a more comprehensive pic¬ 
ture. GRACE densities were taken from a higher altitude 
than the Solomon et al. study. 

The differing results for TIEGCM compared to these 
spacecraft observations may be related the anomalous con¬ 
stituent concentrations found at varying thermospheric 
heights during this deep solar minimum period. Previous 
work has found enhanced concentrations of helium, particu¬ 
larly in the winter hemisphere at ~476km [see Thayer et al., 
2012, and references therein], which is within the range 
of GRACE altitudes during the 2009 period studied here. 
Comparing CHAMP and GRACE observations in 2007 and 
2008, Thayer et al. found that altitude and latitude response 
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Figure 6. Orbit-averaged DTM and TIEGCM with CHAMP DA densities compared with GRACE 
observations for the period beginning 2009 March 15 in the upper row, and 2003 October 15 in lower 
middle row. Density difference between GRACE and the two model runs are shown for the 2009 period 
in the upper middle row, and for the 2003 period in the lower row. The altitude of the GRACE space¬ 
craft during the 2009 period varied between ~506-454 km, and during the 2003 period varied between 
~522-466 km. 


in thermospheric mass density was influenced by the rela¬ 
tive amount of helium (He) and oxygen (O) present. Liu 
et al. [2014] compared CHAMP and GRACE data to the 
semi-empirical NRLMISISE-00 model during the same 2008 
period, finding wintertime helium concentrations exceeding 
model results by 30-70%. 

Liu et al. attributes the high concentrations to the ex¬ 
tremely low He/O transition altitudes in the winter hemi¬ 
sphere in this cold, contracted thermosphere. This is a par¬ 
ticularly important aspect to consider for forecasting pur¬ 
poses since the presence of He instead of O at the same tem¬ 
perature would act to increase the drag coefficient [Thayer 
et al, 2012]. In this work DTM represents this differing 
concentration at GRACE altitudes better than the physical 
model, as is clear from Table 1. DTM is trained on large ob¬ 
servational data sets including the 2007-2008 period, which 
is more representative of the deep minimum conditions than 
data used for training in earlier versions of DTM and the 
NRLMISISE-00 model. The assimilation scheme used with 


TIEGCM only has a local impact at the lower CHAMP al¬ 
titudes during the period being run. 

It is clear from the results that while the semi-empirical 
model runs well during solar minimum, the physical model is 
more accurate when more stormy conditions are introduced. 
TIEGCM generally performs better during solar maximum 
conditions (Figure 2), however DTM also mirrors changes 
in CHAMP densities well. Both models underestimate the 
observed density during the period beginning 2003 March. 
This is similar to the period of increased activity beginning 
2003 October, where underestimations also exist (although 
TIEGCM densities match observations very well, as shown 
in Figures 4 and 6). Sutton et al. [2005] examines the geo¬ 
magnetic impact of the Hallowe’en storms, and also finds the 
NRLMSISE-00 model to underestimate densities in compar¬ 
ison with CHAMP data during times of maximum geomag¬ 
netic activity. 

Since the Hallowe’en storms are such a well-studied event, 
our results from this period can be compared to previous 
studies. Sutton et al. [2005] find density measurements ex¬ 
hibit enhancements of 200-300%. Sutton et al. [2006] uses 
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GRACE and CHAMP data to observe the thermosphere 
neutral density response to the X-class flares on 2003 Octo¬ 
ber 28 and November 4. Sutton et al. finds an increase of 
^50-60% in thermospheric density associated with the first 
flare at low to mid-latitudes, and a ~35-45% increase in re¬ 
sponse to the second flare. Figure 1 of Qian and Solomon 
[2012] shows the density response to the two flares using 
TIEGCM simulations as well as GHAMP observations. A 
^100-200% density enhancement occurs, with storm re¬ 
sponse largest at high latitudes. These results compare well 
with the enhancements shown in Figure 4, with increases of 
up to ~200% observed, and a slightly larger response found 
for the first flare compared to the second. 

It is interesting to note there is a smaller enhancement ob¬ 
served in GRAGE and interpolated model densities in Fig¬ 
ure 6, the observations obtained at a higher altitude (by 
^100 km) than the CHAMP observations. The storm-time 
responses examined here may also be affected by constituent 
composition varying with altitude, as mentioned previously 
for the solar minimum period studied. Liu et al. [2014] found 
the He enhancement causes differing perturbations of mass 
density with altitude after geomagnetic activity. In fact, 
these mass density perturbations are less enhanced in the 
winter hemisphere at GRACE altitudes (^25%) than the 
lower CHAMP altitudes (^60%) during the 2008 minimum. 
This could account for the smaller enhancements observed 
at GRACE altitudes here, however a more thorough inves¬ 
tigation of latitudinal variation would be needed for confir¬ 
mation. 

Although the physical model clearly outperforms the 
semi-empirical model during stormy conditions, a number 
of improvements could still be made to the TIEGCM set-up 
used here. This work focuses on 2-month periods around 
solar equinox, and future work should also examine the sol¬ 
stice periods. However, for longer periods of study it must 
be considered that thermospheric neutral density and com¬ 
position exhibits strong seasonal variation. Here the inter¬ 
action with the lower atmosphere becomes important, how¬ 
ever the standard version of TIEGCM used in this work only 
includes a simplified representation of the lower boundary 
(see Section 2). Qian et al. [2009] found that TIEGCM 
produces a better representation of the annual/semi-annual 
variation with gravity wave parameterisation and an eddy 
diffusion coefflcient included at the lower boundary. The 
standard version of TIEGCM used in this work specifies an 
eddy diffusion coefhcient as a constant at the lower bound¬ 
ary, which decreases exponentially with increasing pressure 
levels. Qian et al. modified this globally uniform value to 
account for seasonal variation, representing the coefficient as 
a Fourier series with four harmonics per year. Research fol¬ 
lowing on from the ATMOP project will include forcing the 
lower boundary of TIEGCM with a non-hydrostatic model 
of the lower atmosphere, which includes gravity wave pa¬ 
rameterisation and non-migrating tides. 

This work has presented validation periods at solar max¬ 
imum and a particularly deep solar minimum. In order to 
gain a complete picture of how solar variability may be af¬ 
fecting the performance of the two models it may be useful 
to examine less deep solar minimum periods, as well as ‘so¬ 
lar medium’ periods that lie midway through the ascending 
and descending phases of a solar cycle. Thermospheric mod¬ 
els are known to be strongly influenced by external drivers, 
including geomagnetic and solar forcing besides boundary 
forcing. Much previous work has found that, alongside other 
improvements to the model, simply adjusting the F10.7 so¬ 
lar flux can improve results [Thayer et al., 2012; Matsuo 
et al, 2013; Liu et al, 2014]. Future work incorporating 
both lower boundary improvements to TIEGCM and exam¬ 
ining results over longer time periods will better highlight 
the possible impact that seasonal and solar variability may 
be having on the model densities. 

It is clear from Figure 6 that interpolating the CHAMP- 
assimilated TIEGCM density results to higher satellite 


height did not have a major impact, with similar errors 
found to previous results. However, with such limited obser¬ 
vations used during assimilation (one data point rather than 
a whole 2D map at each timestep) this is likely highlighting 
the accuracy of the model rather than saying too much about 
the assimilation procedure. As mentioned previously, the 
satellites are at slightly different local times, and enhanced 
He concentrations are found at the higher GRACE altitude 
particularly during the deep solar minimum in 2009. This 
likely introduces additional density errors that may coun¬ 
teract the positive impact of assimilating CHAMP data at 
these altitudes. At the lower altitudes, a positive but lim¬ 
ited improvement has been obtained using data assimilation 
with a general circulation model, with an overall improve¬ 
ment of ~4% found. Including more relevant observations 
in the assimilation procedure would likely increase this per¬ 
centage improvement. 

Future work will also include improvements to the data 
assimilation techniques developed in the ATMOP project, 
particularly implementing incremental analysis updates (see 
Henley et al., manuscript in preparation, 2015, for further 
discussion). The EnOI method developed could be con¬ 
verted to an EnKF system without much difficulty. An EnOI 
system is similar to EnKF, however information from the 
observations which improve the main model does not feed 
back into the ensemble. Assessment of the other versions of 
DTM also developed as part of the ATMOP project could be 
undertaken using the same method as in this paper. DTM- 
2013 is supplemented by 2.5 years of GOCE satellite data, 
uses the 30 cm radio flux as a solar proxy, and 3-hour Am in¬ 
dex as geomagnetic proxy. To improve predictions to 3-days 
out, another method (DTM-nrt) was developed to predict 
temperature corrections to DTM based on a neural network. 
Recent research highlights that the DTM-nrt version of the 
model is more accurate than the DTM-2012 studied here 
for 24-hour forecasts [Choury et al., 2013]. Comparisons 
to and assimilation of other types of observations would 
also be beneficial, for example Fabry-Perot interferometer 
wind measurements, or total electron content. An increase 
in the availability of more real-time in-situ observations of 
the thermosphere will aid this work and space weather fore¬ 
casting in general. Data from missions such as SWARM 
[Frits-Christensen et al, 2006] and Drag and Atmospheric 
Neutral Density Explorer [Pilinski, 2008] may prove useful 
for this purpose. 

With further improvements, the use of general circulation 
models in operational forecasting, in addition to empirical 
models currently used, is certainly plausible. Future work 
will allow near-real-time assimilation of thermospheric data 
into TIEGCM for forecasting. 
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